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Freezing  of  Soil  with  an  Unfrozen  Water  Content 
and  Variable  Thermal  Properties 

VIRGIL  J.  LUNARDINI 

INTRODUCTION 

The  mathematical  theory  of  conductive  heat  transfer  with  solidification  has  been  largely 
confined  to  materials  that  change  phase  at  a  single  temperature.  The  best-known  problem  of 
this  type  is  that  of  Neumann,  and  its  solution  has  been  widely  used  for  the  freezing  of  soils 
(Neumann  1860,  Berggren  1943,  Carslaw  and  Jaeger  1959).  However,  for  media  such  as 
soils,  the  phase  change  can  occur  over  a  range  of  temperatures  (Anderson  and  Tice  1973, 
Tice  et  al.  1978,  Lunardini  1981a).  In  other  words,  at  any  temperature  below  the  normal 
freezing  point,  there  will  be  an  equilibrium  state  of  unfrozen  water,  ice  and  soil  solids.  Figure 
1  shows  the  geometry  for  a  semi-infinite  soil  mass,  initially  at  a  temperature  above  freezing, 
that  freezes  due  to  a  constant  surface  temperature  held  below  the  freezing  point.  The  phase 
change  is  assumed  to  occur  within  the  temperature  limits  of  Tm  and  7f,  representing  mini¬ 
mum  and  maximum  phase  change  temperatures. 

Figure  2  shows  the  unfrozen  water  $  as  a  function  of  temperature  for  a  typical  soil.  At  T{ 
all  of  the  water  is  in  the  liquid  form,  while  at  Tm  the  free  water  is  all  frozen.  There  may  be  a 
residual  amount  of  bound  water,  denoted  by  £f,  that  will  remain  unfrozen  even  at  very  low 
temperatures.  It  will  be  assumed  that  for  T<  Tm,  unfrozen  water  may  exist  but  no  phase 
change  will  occur.  The  region  Tm  <  T  <  Tf  is  called  the  zone  of  phase  change,  or  the  mushy 
zone.  In  this  region,  water  will  solidify  to  ice,  and  unfrozen  water  and  ice  will  coexist.  As 
(Tf  -  T,,,)  —  0,  the  phase  change  will  approach  the  typical  Neumann-type  problem,  which  is 


Figure  1.  Geometry  for  solidification  with  Figure  2.  Unfrozen  water  vs  temperature, 
a  phase  change  zone. 


particularly  valid  for  coarse  materials  such  as  sands  and  gravels.  The  form  of  the  £  function 
for  soils  can  be  expressed  by  different  functional  relations.  The  simplest  relation  is  linear: 


I 
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'l 
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£  =  &  +  -£r  rf)  o> 

where  A£  =  £0-£f  and  A Tm  =  Tf-Tm.  Another  functional  relation,  which  can  closely 
model  the  data  and  is  easy  to  manipulate  analytically,  is  a  quadratic  form: 

£  =  +  (2) 

If  £o,  £f  and  A7-,,,  are  the  same  for  these  functions,  then  the  mean  unfrozen  water  slope  d$/dT 
will  be  identical. 

The  thermal  conductivity  and  the  specific  heat  within  the  mushy  zone  are  functions  of  the 
unfrozen  water  and  may  be  represented  by 

(kf- Au) 

k  =  ku - (3) 

(Cf-Cu) 

c  =  Cu-  (£-£.)  (4) 

where  kf,  ku  =  fully  frozen  and  fully  thawed  thermal  conductivities  and 
Cf,  Cu  =  fully  frozen  and  fully  thawed  specific  heats. 

Obviously  these  properties  are  functions  of  the  particular  form  of  the  unfrozen  water  func¬ 
tion  (Frivick  1980).  Within  the  fully  frozen  region  (zone  1)  it  is  assumed  that  the  thermal 
properties  are  constant  and  equal  to  the  frozen  values,  while  for  the  thawed  region  (zone  3) 
the  properties  are  constant  and  equal  to  the  thawed  soil  values. 

Tien  and  Geiger  (1967)  and  Ozisik  and  Uzzell  (1979)  used  an  unfrozen  liquid  content  that 
varied  with  position  within  the  two-phase  zone.  They  did  not  deal  with  soil  systems,  how¬ 
ever.  Cho  and  Sunderland  (1974)  found  a  solution  for  the  freezing  of  a  material  with  the 
thermal  conductivity  varying  linearly  with  temperature;  however,  the  material  changed 
phase  at  a  single  temperature. 


BASIC  EQUATIONS 

Consider  a  small  volume  of  material  within  the  mushy  zone.  Energy  will  be  conducted  in 
and  out  of  the  volume,  and  latent  heat  will  be  released  during  solidification  (Fig.  3).  Thus  the 
problem  is  one  of  conduction  with  a  distributed  energy  source.  The  governing  equations 
were  derived  by  Lunardini  (1985,  1988).  The  net  conduction  is 


«»-«»+** =  fx(kd^htAx- 


The  latent  energy  released  due  to  solidification  of  a  mass  of  water  A mw  is  as  follows 
A<7g  =  -  f  Amw  =  -( 7dA£Ax. 
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Figure  3.  Heat  flow  in  the 
mushy  zone. 


The  energy  equation  in  the  mushy  zone  is 


The  most  general  case  will  be  a  problem  with  three  regions  as  shown  in  Figure  1 .  The  equa 
tions  for  region  1  are 


d2T,  j_ar, 

dx1  a ,  dt 

U 0.  0  =  Ts 

UXut)  =  Tm 

dT,(Xf  ?M) 

*'  -~d~x  k(X,)  ~W~ 

For  region  2  (the  mushy  zone) 

T,(X„t)  =  Tm 
Ti(X,  t)  =  Tf 

dT,(X)  dT,(X) 

k{X)  ~ k>  ~nr~  • 

The  thawed  zone  is  governed  by 


• 

d}T,  1  dT, 

K 

r. 

dx 1  aj  dt 

f: 

r 

lim  T,  (x,  t)  =  To 

K 

X—  00 

* 

T,(x,o)  =  To 

l 

T,(X,t)  =  T{. 

(8 


(8a; 

(8b; 


(8c; 


(9; 

(9a; 

(9b; 


(9c; 


do; 

(10a 

(10b 


(10c 


These  equations  can  be  partially  nondimensionalized  by  using  the  following  relations: 


T(~T\ 

i  =  1,2,3 

0-  ^ 

0o  — 

To -7} 

T(-Tm 

0  TV-  7-m 

Tf-T^ 

C*J2  <t> 

Sj 

C,o  0 

(Jo  —  « 

ox 

X  —  00  k\2 

Xo  = 

00  ^30 

p 

o 

ii 

c  C,(Tf-Ts) 

T  7d(A£ 

Cfu- 

With  these  relations  the  thermal  conductivity  and  specific  heat  in  region  2,  for  the  linear  £ 
case,  are 


A"  =  *.(1  +0,6:) 

(11a) 

C  =  Cod +0:6:). 

(lib) 

For  the  mushy  zone  with  variable  thermal  properties,  eq  1  la  and  b  will  use  k„  =  ku  and  C0  = 
Cu.  However,  the  use  of  the  general  notation  k0  and  Cu  allows  the  properties  in  region  2  to  be 
specified  independently  of  the  frozen  or  thawed  values  if  the  thermal  properties  are  constant 
(Appendix  A). 

In  the  mushy  zone  the  following  transformation  will  be  used: 

l  6‘ 

0  =  7H  k(6{)d6L 

Ko  0 

(12) 

For  the  linear  £  case  the  function  0  can  be  evaluated  explicitly: 

0  =  d:  +  ^  el 

(12a) 

Then 

k  =  ko  Vl +2/3,0 

(13) 

„  VT+20T0-1 

02  =  0, 

(14) 

The  transformed  equations  are  as  follows: 

d%  J_  dfl. , 

dx !  ai  dt 


U5) 


^,(0,0  =  0 


(15a) 


*.  yx  {x"l)  =  *<*■>  5?  (*»')• 


(15c) 


The  equations  for  region  2  are  written  only  for  the  linear  £  assumption  (the  details  are  given 
in  Appendix  A): 

-  a 

(16) 

UX,t)  =  o 

(16a) 

i(Xut)  =  1+/3./2 

(16b) 

ko  (X,t)  =  k,  ^  (X,t) 

(16c) 

«0 F  =  (1  +  <*-&,)*  +  (1  +  2/3, *)v\ 

(17) 

The  thawed  region  has  the  following  relations: 

320,  _  1  30, 
dx1  a>  dt 

(18) 

lim  6,(x,t)  =  -0o 

X—  oo 

(18a) 

0,(*,O)  =-</>„ 

(18b) 

© 

II 

•W* 

* 

(18c) 

Exact  solutions  for  the  three-zone  problem  with  variable  thermal  properties — even  with  a 
linear  £  function — have  not  been  found.  It  is  possible  to  obtain  approximate  solutions  using 
the  heat  balance  integral;  however,  before  doing  this  it  is  useful  to  examine  the  simpler  two- 
zone  problem.  If  the  surface  temperature  Ts  is  greater  than  or  equal  to  the  minimum  phase 
change  temperature,  then  a  completely  frozen  zone  will  not  exist.  Thus  we  need  only  examine 
regions  2  and  3. 


TWO-ZONE  PROBLEMS 

The  two-zone  problem  is  simpler  than  the  three-zone  case  and  will  lead  to  results  that  can 
simplify  the  need  for  the  full  three-zone  problem.  The  linear  unfrozen  water  case  will  be  ex¬ 
amined  for  both  variable  and  constant  thermal  properties,  while  the  quadratic  water  content 
case  will  be  evaluated  only  for  the  constant  property  problem.  This  will  be  shown  to  be  ade¬ 
quate  for  the  general  problem. 

Linear  unfrozen  water  function 

Variable  thermal  properties 

Equations  16-18  are  valid  for  this  case  except  that  the  boundary  condition  in  eq  16b  be¬ 
comes 


5 


to  vtHr.xr  vw  kwk\ 


i(0,0  =  4>+ ^  =  P.  (18d) 

An  approximation  of  the  solution  may  be  obtained  with  the  heat  balance  integral  method. 
The  heat  balance  integral  method  has  been  adapted  to  problems  of  freezing  in  soils  systems 
by  Lunardini  (1981b,  1982,  1983)  and  Lunardini  and  Varotta  (1981).  The  equations  for  the 
heat  balance  integral  method  are  well  known  and  will  not  be  derived  here;  the  interested 
reader  can  consult  Lunardini  (1981a)  for  details.  Referring  to  Figure  1  and  Appendix  B,  eq 
16  and  18  become 
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X  X 

\  VTTW  ^  dx  =  ±[\  F(x,t)dx-F>X) 

i[\ 


For  the  heat  balance  integral  treatment,  eq  1 8a  and  b  become 


0a(8,O  =  ~<t>  o 


(6,0  =  0. 


Quadratic  temperature  profiles  are  assumed  for  ^  and  6,  since  experience  has  shown  that 
they  yield  good  results  for  the  heat  balance  integral  method: 

t  =  b(X-x)  +  c(X-xY  (21) 


cX 2  =  P-bX. 


To  simplify  the  algebra  the  following  parameters  are  defined: 


X  =  2  y\!a}t 


5-X  =  BX. 


The  use  of  eq  23  and  24  is  not  a  constraint  on  the  final  solution. 


?! 


The  solution  of  eq  19  and  20  is  straightforward  but  tedious.  The  algebraic  manipulations 
are  given  in  Appendix  B.  The  unknown  parameters  7  and  B  can  be  found  from  the  following 
equations: 


«(■?  +  ') 

KQ ,  fl(y  +  l)  =  [f,  (y  +  y)  +  F(Q2  -  1) 


(25) 


(26) 


JlfaK  +  0>(2K  +  A) 


sfN  A(N- 1) 


Q,  =  , _ 

\[20JC  +  0,A 

N  =  1  +  2/3, P 


.  2Xo 
~B 


K  =  P-/1 


F,  =  —  (I  +  <7.  -  ft,), 
ao 

Constant  thermal  properties 

The  constant-thermal-properties  solution  follows  from  the  preceding  case  if  /32,  —0  and 
0,  —  0.  It  can  be  shown  that,  for  /32,  =0, 


lim 

0i  —0 

Q,  =  1 

(27a) 

lim 

0,-0 

fit  =  1 

(27b) 

lim 

0i  —0 

i  =  0j. 

(27c) 

For  the  constant-property  case,  a0  becomes  a,  and  =  a,  Xo  =  X.  Then  the  solution  is 
(p  -  j)(B  +  1)  -  (l  +  ^)(1  +  <7)a„  =  0.  (28) 

The  parameter  7  is  again  found  from  eq  25. 


7 


For  this  case  an  exact  solution  is  possible  by  using  a  similarity  transformation,  as  was 
shown  by  Lunardini  (1985).  The  solution  is 


h. 

<t> 


erf/x/2 

\  VI  +  of 

erfyVa}j(l  +  a) 


h 


erfcy 


XerfyVa»(l+u)  =  VaJ2(l  +  a) erfcy  e-'^0'”*1  +  o]~ 11 


(29) 


(30) 

(31) 


where 

,  * 

erfx  =  — :  J  c~z’  dz 

vx  0 

is  the  familiar  error  function. 

Quadratic  unfrozen  water  function 

When  a  quadratic  unfrozen  water  function  is  used  (for  example,  eq  2),  it  is  not  possible  to 
find  a  closed-form  solution  for  variable  thermal  properties.  Thus  a  heat  balance  approxima¬ 


tion  will  be  used  for  constant  thermal  properties. 

The  equations  for  region  2,  using  eq  2,  are 

«a  IF  =  Tt  N1+2<^-^n  <32> 

6,(X,t)  =  0  (32a) 

0a  (*„0  =  1  (32b) 

kl  ^x(X,t)  ~  ki  lix  (X,t)‘  (32c) 

The  heat  balance  integral  form  of  eq  32  is 

x 

«a  (X,0-  Yx  <°*'>)  =  ^  f  10  +  2<r)0a  -  o&l]dx.  (33) 

Equation  20  is  still  valid  for  region  3. 


The  temperature  profile  for  region  3  is  again  eq  22.  The  temperature  in  region  2  is  assumed 
to  be 


02  =  b,(X-x)  +  c,(X-xY 


(34) 


where 


.•  A'V.VJtVA' 


-  *-  |  . 

The  equation  for  B  is  now 

(<t>B  -  2\)(B  +  3)  =  +  <t>)(l  +2a-  -  J-  -?£-] .  (35) 

The  parameter  y  is  again  found  from  eq  25. 

The  two-zone  solutions  can  be  compared  by  considering  some  specific  cases.  For  example, 
consider  a  typical  soil  with  the  following  properties  suggested  by  Nakano  and  Brown  (1971): 

T„  =  4°C 

Tt  =  -4 

T,  =  0 

Tm  =  -4  (thus  0  =  0o  =  1) 

£o  =0.20gWa^ 
gram  solids 

if  =  0.0782 

7d  =  1.68  gram  solids/cm 3 
l  =  80  cal/gram  water 
ST  =  0.1539. 

The  soil  thermal  properties  and  the  results  of  several  cases  are  summarized  in  Tables  1  and  2. 
Cases  1-3  show  that  the  effect  of  specific  heat  variation  is  not  important  and  can  safely  be 
neglected.  However,  case  4  indicates  that  the  thermal  conductivity  can  cause  15-25%  varia¬ 
tions  in  the  rate  of  growth  of  the  freezing  zone.  Case  5  uses  average  values  of  k  and  C  within 

Table  1.  Effect  of  thermal  properties  on  freeze  of  soil  with  average  properties  and  linear  i. 


c. 

(cal/cm'-'C) 


k, 

fcal/s-cm- °C) 


Difference 
from  case  l 
<*> 


Comment 


1  0.63 

0.0058 

0.431 

-0.1429 

0.3988 

— 

Base  case  with  variable  k,  C. 

! 

2  0.63 

0.0058 

0.431 

0 

0.3996 

0.2 

Constant  specific  heat,  C  =  0.63. 

j 

3  0.54 

0.0058 

0.431 

0 

0.4016 

0.7 

Constant  specific  heat,  C  =  0.54. 

i 

4  0.54 

0.0083 

0 

0 

0.4575 

14.7 

Constant  k,  C.‘ 

’ 

5  0.585 

0.0071 

0 

0 

0.4126 

3.5 

Constant  k,  C.  t 

6  0.54 

0.0083 

0 

0 

0.4277 

7.2 

Constant  k,  C,‘  exact  solution. 

■ijrtsi  wi, 


Table  2.  Effect  of  thermal  properties  on  freeze  of  soil  with  extreme  property  variations. 

Difference 


Case 

C. 

(cal/cm‘-°C) 

k, 

(cal/s-cm-  °C) 

B, 

0, 

7 

/rom  case  / 

Comment 

1 

0.63 

0.0058 

1 

-0.50 

0.4626 

Extreme  properties,  linear  £,  variable  k,  C. 

2 

0.63 

0.0058 

1 

0 

0.4607 

-0.4 

Constant  specific  heat,  C  =  0.63,  linear  {. 

3 

0.315 

0.0058 

1 

0 

0.4696 

1.5 

Constant  specific  heat,  C  =  0.315,  linear  { 

4 

0.315 

0.0116 

0 

0 

0.5671 

22.6 

Constant  k,  C,*  linear  £. 

5 

0.473 

0.0087 

0 

0 

0.4726 

2.2 

Constant  At,  C,  t  linear  {. 

6 

0.315 

0.0116 

0 

0 

0.5304 

14.7 

Constant  k,  C,  •  linear  £,  exact  solution. 

7 

0.315 

0.0116 

0 

0 

0.5226 

13.0 

Constant  4r,  C*.  quadratic  (. 

1 


KrCn 

BpS  l 


* 


e.  <t>  =  0.1,  an  =  0.1,  0.3. 


f.  <t>  —  0.1,  fn  —  1.0,  0.6. 


Figure  4  (cont’d). 


THREE-ZONE  PROBLEMS 

Since  the  variable-property  case  can  be  adequately  handled  by  an  appropriate  constant- 
property  solution,  only  the  constant-property  problem  will  be  examined. 

Linear  unfrozen  water  function 

Exact  solution 

Equations  15  and  18  are  the  governing  equations  for  regions  1  and  3,  while  the  equations 
for  region  2  are 


d 


02(X,t)  =  0 


Oi(X„t)  =  1 


The  similarity  method  was  used  by  Lunardini  (1985)  to  obtain  an  exact  solution  to  this 
problem.  The  solution  is  as  follows: 


=  (1  -  0)11  + 


\2\Jait 


erf7v«J2(l  +  o) 


-  erf (—4=J\ 


erf7Vo;32(l  +  a) -erfr?vai2(l  -4-  cr) 


(2X7/) 

.  c  _  * 


erfc  7 


X,  =  2rjVa,7  (40) 

X  =  27VS7 T.  (41) 

The  parameters  ij  and  7  are  found  from  the  simultaneous  solution  of  the  two  following  equa- 


Vajj(l  +  a)  (erfc7)e^'(0,’■(,  +  =  x[erf7V^7(TT7)  -  erfWa„(l  +  a)]  (42) 

(0-  De'”'!1  -“'■(l+<’)l[erf7V7ri(l  +<t)  -  erfi,Va,2(l +ff)]  =  *„Va„(l +a)erfij.  (43) 

Lunardini  (1985)  showed  that  this  solution  approached  the  Neumann  solution  as  the  phase 
change  zone  decreased  (Table  3).  It  is  clear  from  Table  3  that  the  solution  does  converge  to 
the  Neumann  solution  as  ( Tf -  Tm)  —  0.  The  thaw/freeze  interface  greatly  exceeds  the  value 
for  the  Neumann  solution  if  phase  change  occurs  over  a  4°C  zone  as  in  the  example  calcula¬ 
tion. 

Heat  balance  integral  solution 

The  heat  balance  integral  equations  for  the  three-zone  problem  are  as  follows: 


■IS  <*.«>  - 


MO,  t)  =  <t> 


i  «•'>  ■  »• 


Quadratic  temperature  profiles  for  the  three  regions  lead  to 

*■  *  *[(J3?)'  - 


9,  =  1  +  bI(Xl-x)  +  c1(Xl-x)1 
&i  =  e,(X-x)  +UX-xY 


where 

b>X  =  2k 


a-i) 

CiX\[  -Rf  =  0  -  1  -  2*„(1  ~R)( 4  -  4) 


(46c) 


(47) 


(48) 

(49) 


< 

S 

5 


I 

p- 


s“ 

v' 


V 


* 


« 


1 

-4 

0.0617 

0.3029 

33.33 

8.13 

25.2 

2 

-2 

0.1135 

0.2576 

28.34 

14.95 

13.39 

3 

-1 

0.1376 

0.2272 

25.0 

18.12 

6.88 

4 

-0.5 

0.1492 

0.2106 

23.17 

19.65 

3.52 

5 

-0.1 

0.1571 

0.1946 

21.41 

20.69 

0.72 

Neumann 

0 

0.1606 

— 

21.15 

21.15 

0 

*  For  /  =  24  hours. 

T.  =  4<>C  Tt=  -6  rf  =  0°C  k,  =  0.00828  cal/s-cm-°C 
*,  =  0.00703  k ,  =  0.00578  C,  =  C,  =  C.  =  0.165  cal/cm>-“C 
^  =  0.0605. 


The  solutions  to  eq  44-46  are  straightforward;  the  details  will  be  omitted.  The  results  are 
given  below: 


=  (1  -RYy'a, 

(*-  1)(3-t,2)-(6  +  >j2)|(1  ~R)k} 
2\R  (1  +  a) 


1  - 


B 


'(«  ~  s)l  '  0 

Rany'(l-2R+  =  0. 


(50) 


(51) 

(52) 

(53) 


These  four  equations  can  be  easily  solved  for  rj  and  7. 

Table  4  shows  that  the  approximate  solution  is  in  error  by  less  than  6%  when  compared  to 
the  exact  solution. 


Table  4.  Comparison  of  exact  and  heat  balance  integral  solu¬ 
tions  with  linear  $  and  constant  k  and  C. 


Heat  balance 


Case 

Tm 

(°C) 

Exact  solution 

integral  solution 

Variation  (%) 

V 

7 

V 

7 

V 

7 

1 

-4 

0.0617 

0.3029 

0.0622 

0.2938 

0.8 

-3.1 

2 

-2 

0.1135 

0.2576 

0.1150 

0.2411 

1.3 

-6.4 

3 

-1 

0.1376 

0.2272 

0.1390 

0.2168 

1.0 

-4.6 

4 

-0.5 

0.1492 

0.2106 

0.1505 

0.2050 

0.9 

-2.6 

5 

-0.1 

0.1571 

0.1946 

0.1595 

0.1958 

1.5 

0.6 

T.  =  4°C  T,  =  -6  Tf  =  0°C  *,  =  0.00828  cal/s-cm- °C 

k,  =  0.00703  k,  =  O.OOJ78  C,  =  C,  =  C,  =  0.165  cal/cm’-°C 
Sr  =  0.0605. 


Quadratic  unfrozen  water  function 

The  equations  for  the  quadratic  unfrozen  water  relation  are  identical  to  eq  44-46  except 
that  the  energy  equation  for  the  mush  region  (eq  45)  is 

x 

{X' 0  ~  d~t  {X"  °]  =  ft  J  [(1  +  2a)Bl  -  +  (1  +  a)X'  ■  (54) 

The  quadratic  temperature  approximations  (eq  47-49)  are  used  again.  The  solution  is 
again  given  by  eq  50-52,  with  the  mushy  zone  equation  given  by 


1  -  2z  -  R1  ^  ^  (2z’-7z+18)j  =0  (55) 

where 


Case  1  of  Table  4  was  evaluated  for  the  quadratic  £,  and  it  was  found  that  j/  =  0.0572  and 
y  =  0.2730.  These  values  differ  from  the  linear  £  approximation  by  about  8%. 

The  quadratic  £,  three-zone  problem  is  a  function  of  X,  a,  and  ku.  Graphical  so¬ 

lutions  are  shown  in  Figure  5  for  typical  soil  parameters. 


CONCLUSIONS 

The  mathematical  model  used  assumed  the  latent  heat  to  be  a  source  of  energy  distribution 
throughout  the  volume  of  a  soil  with  phase-change  temperature  limits  of  Tf  and  Tm.  This 
contrast  with  the  treatment  of  the  latent  heat  is  totally  released  at  the  upper  phase-change 
temperature  Tf.  A  comparison  of  the  exact  solution  for  the  former  case  showed  that  it  con¬ 
verged  exactly  to  the  Neumann  solution  as  T(  -  Tm  approached  zero.  Thus,  the  mathe¬ 
matical  model  is  based  on  sound  physical  principles. 

The  existence  of  a  mushy  zone  can  have  a  significant  effect  on  the  mechanical  strength  of 
freezing  or  thawing  soils.  For  many  predictions  of  the  strength  of  a  freezing  soil,  the  assump¬ 
tion  is  made  that  the  soil  temperature  is  that  calculated  from  the  Neumann  solution.  A  soil 
with  a  mushy  zone  would  have  a  completely  frozen  layer  of  soil  that  is  thinner  than  that  for 
the  Neumann  case  and  therefore  would  have  less  bearing  capacity.  In  addition  a  thick  zone 


of  frozen  soil  would  exist  that  has  a  variable  unfrozen  water  content.  Again  the  presence  of 
this  liquid  water  will  decrease  the  mechanical  strength  of  the  soil.  The  thawing  case  would 
also  tend  to  yield  a  bearing  strength  of  soil  that  is  less  than  that  predicted  with  the  Neumann 
solution  temperatures. 

The  effect  of  variable  specific  heat  on  the  rate  of  phase  change  is  negligible  for  the  cases 
examined;  thus,  it  is  acceptable  to  use  an  average  specific  heat  value  in  the  mushy  zone.  The 
variation  of  the  thermal  conductivity  with  water  content  is  much  more  significant  and  can 
cause  a  15%  change  in  the  rate  of  freezing  of  the  soil.  However,  it  was  shown  that  the  con¬ 
stant-property  solution,  with  average  values  of  the  thermal  properties  in  the  mushy  zone, 
gives  a  solution  that  is  quite  close  to  that  for  the  variable-property  solution.  It  is  acceptable 
to  use  the  much  simpler,  constant-property  solution  with  average  thermal  properties  to  com¬ 


pensate  for  the  actual  variable  thermal  conductivity. 


Phase  Chonge  Parameter 


5 

*• 


*.  *.  <(. 


A  series  of  graphs  are  presented  of  the  constant-property,  three-zone  problem  for  typical 
ranges  of  soil  parameters.  These  graphs  allow  rapid  predictions  to  be  made  for  the  freezing 
of  soils  with  an  unfrozen  water  content  that  is  a  function  of  temperature  and  variable  ther¬ 
mal  properties. 
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APPENDIX  A:  DERIVATION  OF  THE  MUSHY  ZONE  EQUATION 

The  energy  equation  in  the  mushy  zone  (region  2)  is 


The  dimensionless  temperature  in  region  2  is  defined  as 

„  _  lizh 

1  T(~Tm 

Equation  1  transforms  immediately  to 

a  ( ,  do2\  _  de 2  (y<i  dt  aet 


±  (k  - 
dx  \  dxj 


(y<i  di  ae, 

(Tf-Tm)de>  at  ■ 


With  the  parameters  defined  in  the  text,  it  is  also  clear  that,  for  the  linear  £  case, 

A:  =  Aru(l  -t-  0,02)  (A4a) 

C  =  Cu(l  +  0,Ot).  (A4b) 

Equations  A4a  and  b  are  written  explicitly  for  the  definitions  of  k  and  C  discussed  earlier. 
However,  it  will  be  advantageous  to  allow  for  the  case  where  the  mushy  zone  will  have  con¬ 
stant  thermal  properties  different  than  those  of  the  thawed  state.  Thus,  we  will  define  k0  and 
Co  to  be  any  constant  values  of  the  thermal  conductivity  and  specific  heat  desired.  Then 

k  =  Ar0(l  +  /3,02)  (A5a) 

C  =  C0(l  +  0202).  (A5b) 

Clearly,  for  the  mushy  zone  with  variable  thermal  properties,  k0  =  ku  and  C0  =  Cu.  How¬ 
ever,  if  we  wish  to  examine  the  constant  thermal  property  case,  then  we  can  simply  set  k0,  C0 
to  any  convenient  values.  For  example  we  could  let  the  properties  of  the  mushy  zone  be  the 
mean  values  of  the  frozen  and  thawed  states  and  define  k0  and  C0  appropriately. 

If  the  thermal  conductivity  varies  with  temperature,  the  Kirchoff  transformation  can  be 
used  to  define  a  new  temperature  variable: 

1  9‘ 

i  =  Y  j  k(y)dy  (A6) 

Ko  o 

where  y  is  a  dummy  variable. 

From  eq  A6 


The  derivation  will  proceed  for  the  case  of  a  linear  £  function.  The  quadratic  assumption  can 
also  be  used,  of  course,  but  the  equations  will  be  more  complicated. 


From  eq  A4  and  A6  an  explicit  relation  between  and  0j  can  be  found: 

4>  =  ei+  y  02.  (A9) 


It  follows  that 


k  =  AfoVl  ■+■  2(3\\p 
C  =  Co(l  -  ft,  +  ft.x/l  +2fttf). 


Equation  A3  then  becomes 


fc  vrraj -gS- -  (i 

Equation  A12  can  be  put  into  a  more  convenient  form  by  noting  that 


(A  10) 
(All) 


(A  12) 


, _  d\b  Id  u 

Vl+2ft*  10+ 2ft*)4].  (A  13) 

Finally  the  energy  equation  for  the  mushy  zone  with  a  linear  unfrozen  water  content  function 
is 


«oVl  +2/3,1^  4^-  =  4  f(l  -  ft,  +  *«)*  +  (1  +  2/3,^)Vi  • 


(A  14) 


*1*1 

ij 


<■  V  ' AA  -Jk  •vy’Vli 'J.'  i  »  V*  .  "  W VYL"»  WWW-V ,  «.  S."  *.  v.'i  V  «.7 K"  R_"  K*X  'in-  »i  . 


APPENDIX  B:  SOLUTION  OF  THE  TWO-ZONE  PROBLEM  WITH 
A  LINEAR  £  AND  VARIABLE  THERMAL  PROPERTIES 

The  equations  for  zones  2  and  3  are 


- - d*\P  dF 

'^  +  2S3^  aX2  ~  d, 


mo  =  4>+  ~  =  p 


4(X,t)  =  0 


*•  H<*'>  ’  *•  S  **•» 


gg,  =  j_  aa, 

dxJ  aj  3/ 

6,(6,  t)  =  -0o 


I'  <*•')  - 0 


where 


F  =  F,yp  +  F2(l  +  2/3,^) Vj 


F,  =  —  (1  +  a,  -  ft,) 

CKo 


02  I 

3ao/3i 


The  heat  balance  integral  method  (HBI)  integrates  the  energy  equation  over  the  volume  of 
interest.  For  region  2 


IVTT2W dx-  J  jjd,. 


Using  Leibniz’s  rule  this  is 


r  , - a  r  CIA 

f  VT+20^  £  dx  =  —  j  F\x,t)dx  -  F(X,t)  -  . 


F(X,t)  =  F1 


(B6) 


and  the  HBI  equation  for  region  2  is 


x  x 

\J\+2Md£dx=  ft\\F(x,t)dx-F2x\. 


The  HBI  equation  for  region  3  is  derived  in  the  same  manner  and  is 


-“■5  «*•<>- 


Quadratic  temperature  profiles  for  regions  2  and  3  that  satisfy  the  boundary  conditions  are 
t  =  b(X-x)  +  c(X-xy  (B9) 


—  4>  o 


where 


cX1  =  P-  bX. 


The  moving  interfaces  X  and  5  are  assumed  to  have  the  forms 


X  -  2y\ffXit 


5  =  (B+  Q* 


where  B  is  a  constant.  Using  eq  BIO  with  eq  B8  leads  to  the  following  equation: 


£x  = 


The  solution  to  eq  B13  follows  easily  with  the  aid  of  eq  11  and  12: 


F  /D  \  •  V 

Equation  B7  can  be  solved  to  yield  an  algebraic  equation  for  the  unknown  quantity  B.  From 
eq  B9  we  note 


where  c  is  only  a  function  of  time.  Then  eq  B7  is 


x  lx 

■I  Vi +2 Mdx  =  4r\\nx, 


t)dx  -  FiX\ 


A  4 


■rjrjr. 


'vw - %.<\.  -v^.1- « \  w f.T.v.v.w.v.’.vjv  ,jtxm,wv 


The  solution  of  eq  B16  is  straightforward  but  rather  tedious.  The  results  are 
X 

f  Vi  +  2iWdx  =  Q,X 

o 

x  x 

j  Fdx  =  F,  (y  +  y)  A"  +  F,  j  (1  +  2^)/ldx 


f  (1  +  2t3,t)v‘dx  =  Q2* 


where 


Cl,  v«  + 


\j2fiiK  +  (h(2K  +  A) 
sf2KK+ptA 


3  1-dVlY 

♦  W- >* 


F  =  P  -  /t. 


The  equation  for  fl  is  then 


KQtB  (y  +  l)  =  [f,  +  yj+  Fi(Qi  -  1)  a,. 


A  facsimile  catalog  card  in  Library  of  Congress  MARC  format  is  repro¬ 
duced  below. 
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